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1.  Introduction 


In  studies  of  the  earth’s  anomalous  j^ravity  field,  empirical  covariance 
functions  are  of  great  interest.  The  behavior  of  the  gravity  field  is  reflected 
in  these  functions.  The  magnitude  of  the  variations  and  the  roughness  of  the 
field  is  described.  This  kind  of  information  is  important  and  has  to  be  taken 
into  account  when  gravity  field  related  quantities  are  estimated  from  a  set  of 
observations.  The  method  of  least  squares  collocation  (Moritz,  1980),  is  widely 
used  for  this  purpose.  The  method  needs  a  covariance  function  to  express 
the  relations  between  the  observations  and  between  the  estimated  quantities 
and  the  observations.  The  best  least  squares  agreements  with  the  true 
potential  is  obtained  when  the  empirical  covariance  function  is  used. 

Global  empirical  covariance  values  were  estimated  from  gravity 
observations  by  Tacherning  and  Rapp  (1974).  They  derived  some  covariance 
function  models,  and  an  expression  for  the  global  empirical  covariance  function 
was  found. 

When  studies  of  the  gravity  field  takes  place  in  local  areas,  the  use  of 
high  degree  spherical  harmonic  approximations  is  very  important.  Estimations 
of  gravity  field  related  quantities  are  carried  out  relative  to  the  spherical 
harmonic  approximations  using  the  residual  observations  and  the  local 
empirical  covariance  function.  This  procedure  corresponds  to  the  stepwise 
collocation  (Moritz,  1980),  where  the  solution  of  the  first  step  is  known. 

The  determination  of  a  local  empirical  covariance  function  was  discusseci 
by  Goad  et  al.  (1984).  They  arrived  at  the  following  definition  of  a  local 
covariance  function:  "A  local  covariance  function  is  a  special  case  of  a  global 
covariance  function  where  the  information  content  of  wavelengths  longer  than 
the  extent  of  the  local  area  has  been  removed,  and  the  information  outside, 
but  nearby,  the  area  is  assumed  to  vary  in  a  manner  similar  to  the 
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information  within  the  area".  They  suggest  that  the  Tscherning/Rapp  models 
are  used  in  order  to  adjust  an  expression  to  empirical  values  in  a  similar  way 
as  in  the  global  case.  This  was  done  by  Knudsen  (1987)  using  an  iterative 
least  square  inversion  procedure.  The  empirical  covariance  values  were 
determined  from  altimeter  and  gravity  data  relative  to  a  spherical  harmonic 
approximation  of  degree  and  order  180. 

A  different  approach  is  to  use  the  known  part  of  the  anomalous  density 
distribution  in  the  earth.  The  effect  from  the  topography  is  a  considerable 
part  of  the  gravity  field  and,  as  far  as  the  topography  is  known,  this  effect 
can  be  used  when  the  gravity  field  is  modelled.  Results  obtained  by  Porsberg 
(1984)  show  that  the  topographic  effect  is  a  large  part  of  the  short 
wavelength  part  of  the  gravity  field.  Furthermore,  the  results  show  that  the 
gravity  field  becomes  much  more  homogeneous  when  terrain  reduced  quantities 
from  different  local  areas  are  evalueated.  In  these  studies,  a  density  contrast 
of  2670  kg/m^  was  assumed,  but  actual  density  contrasts  could  also  be  taken 
into  account.  This  was  done  by  Siinkel  et  al.  (1987)  when  they  computed 
terrain  effects  using  both  a  digital  terrain  model  and  a  digital  density  model. 

Covariance  functions  are  needed  in  estimation  methods  like  least  square 
collocation  and  they  have  to  be  consistent  with  the  observations  that  are 
used.  Unfortunately,  the  estimations  of  such  functions  is  hampered  by  the 
lack  of  data  in  the  areas  where  the  gravity  field  modelling  is  to  take  place. 
Consequently,  a  covariance  function  model  has  to  be  chosen.  It  is,  therefore, 
important  to  understand  the  variability  of  the  local  empirical  covariance 
functions.  The  results  mentioned  above  indicate  that  this  variability  decreases 
remarkably  when  the  terrain  effect  is  taken  into  acount.  The  residual  gravity 
field  becomes  more  homogeneous  and  the  choice  of  a  covariance  function  model 
becomes  easier. 
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[n  this  report,  techniques  for  the  estimation  of  empirical  covariance 
functions  are  described.  In  order  to  show  the  importance  of  terrain 
reductions,  previous  results  are  summarized.  Problems  with  centering 
observations  are  treated  and  a  method  for  the  determination  of  a  covariance 
function  model  is  described.  New  studies  of  the  variability  of  the  local 
empirical  covariance  functions  are  carried  out  in  ocean  areas  using  altimeter 
data  along  with  bathymetry  data.  Terrain  effects  are  calculated  from  the 
bathymetry  and  the  role  of  an  isostatic  compensation  of  the  ocean  bottom 
topography  is  evaluated. 

2.  The  Local  Empirical  Covariance  Function 

In  this  chapter  the  estimation  and  modelling  of  a  local  empirical  covariance 
function  are  evaluated.  The  importance  of  a  complete  removal  of  wavelengths 
longer  than  the  extent  of  the  local  area  is  pointed  out,  and  a  method  for 
doing  this  consistently  using  different  kinds  of  observations  is  given. 
Problems  in  using  space  domain  or  frequency  domain  methods  for  the 
estimation  of  a  covariance  function  are  treated,  and  the  effect  of  noise 
associated  with  the  observations  is  evaluated.  Finally,  an  iterative  least 
square  inversion  procedure  for  adjusting  a  covariance  function  model  to  fit 
empirical  covariance  values  is  described. 

2.1  The  local  covariance  function 

In  the  following  AT  is  the  anomalous  potential  T  where  the  information 
content  of  wavelengths  longer  than  the  extent  of  a  local  area  (^i  *^2  )> 

Tref,  13  subtracted.  L  and  L'  are  two  linear  functional  associated  with  the 
observations  y=L(AT)  and  y'=L'(AT)  located  at  (♦,X)  and  (♦',X').  If  the 
averages  of  y  and  y'  over  the  area  are  zero,  then  the  local  autocovariance  of 
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y  (if  L  =  L')  or  the  local  croaacovariance  between  y  and  y'  is  given  by 


♦a  Xj  2Tt 

C(t^)  =  J  j  1  2^  j  yy  **  **  co8(^)  d^dX  (2.1) 

♦  i  X»  o 

where 

coa  V  =  sinf  sint'  +  cos4  cos4'  co8(X-X')  (2.2) 

A  ia  the  size  of  the  area  on  a  unit  sphere, 
a  is  the  axinuth. 

This  rotational  invariant  representation  of  the  covariance  function  is 
calculated  as  an  average  of  the  products  yy'  over  the  area  (the  two  outer 
integrations)  and  an  average  over  the  azimuth  (the  inner  integration). 
(Correaponding  to  homogeneity  and  isotropy  respectively  in  a  plane.)  The 
integration  over  the  local  area  ia  restricted  to  the  observations  y.  When 
goes  from  zero  to  n  the  observations  of  y'  are  used  located  over  the  whole 
sphere.  (See  also  Tscherning,  1985,  and  Sanso,  1986). 

In  practice  the  observations  are  given  in  discrete  points  in  the  area  and 
the  calculation  of  the  covariance  function  is  then  done  by  numerical 
integration  (Tscherning  and  Rapp,  1974).  If  each  observation  y^  represents  a 
small  area  A|  and  yj  represents  an  area  Aj  then 

Cfc  =  I  (2.3) 

L  A^Aj 

With 


fk-i  fij  ^  Vn  (2-4) 

If  the  area  is  subdivided  into  small  cells  holding  one  observation  each  and 
A I  and  Aj  are  assumed  to  be  equal  then  equation  (2.3)  reduces  to 


c.  = 

Nil 


(2.5) 


where  N|,  is  the  number  of  products,  yiyj,  in  the  Ic'th  interval. 

Suppose  T  is  expanded  into  spherical  harmonics,  and  a  global  gravity 
potential  approximation  up  to  degree  and  order  N,  Tn,  is  subtracted  in  order 
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to  obtain  AT.  Then  the  covariance  function*  K{ip),  associated  with  (T-T^)  is 
expressed  by  a  sum  of  a  series  of  Legendre  polynomials  of  order  i*  Pj* 
(Tscherning  and  Rapp,  1974,  and  Tscherning,  1986) 

m)  =  ?  e<(T,T)(-^)‘''‘  P<(cos  1&)  +  I  e,(T,T)(Sir)‘''*  P,(cos 
»=o  rr  i=M+t  “ 

(2.6) 

where 

ei(T,T)  are  the  potential  error  degree-variances, 
ej(T,T)  are  the  potential  degree- variances, 
r,r'  are  the  radial  distances  of  y  and  y' , 

R  is  the  mean  earth  radius,  and 

Rg  is  the  radius  of  a  Bjerhanmar  sphere  (R$  ^  rr'). 

The  integer  N,  relative  to  the  size  of  the  local  area,  is  supposed  to  fulfill  the 
condition  that  27r/N  is  smaller  than  the  extention  of  the  area  where  the  local 
covariance  function  is  determined.  (T-Tm)  is  equal  to  aT  if  an  exact  agreement 
between  T^ef  ®od  T^  exists.  Consequently,  the  error  degree-variances  are 
zero.  The  degree-variances  are  positive  numbers  and  are  related  to  tha 
spherical  power  spectrum  of  the  earth  gravity  field.  It  is  well  known  that  the 
degree-variances  tend  to  zero  somewhat  faster  than  i~’  and  that  the 
Tscherning/Rapp  (1974)  model  4  is  a  reasonable  choice  for  a  degree-variance 
model, 

<r,(T,T)  =  A/((i-l)(i-2)(i+24))  (2.7) 

where  A  is  a  constant  in  units  of  (m/a)*.  For  geoid  heights,  C,  and  gravity 
anomalies,  Ag,  we  then  have  the  degree  variances  associated  with  the 
respective  auto-  and  cross-covariance  functions  (y  is  the  normal  gravity) 

=  l/(r7')e,(T,T)  (2.8) 

<^»(Ag,Ag)  =  (i-l)*/(rr')e,  (T,T)  (2.9) 

<r,(f,Ag)  =  (i-l)/(7r')«f,(T,T)  (2.10) 

See  Tscherning  (1976)  for  degree-variances  associated  with  other  gravity  field 

related  quantities. 
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In  a  flat  earth  approximation  the  sphere  is  replaced  by  a  plane  and 
is  replaced  by  (x,y,s)  where  s‘=x*-»'y*.  The  power  spectral  density,  or,  more 
loosely,  the  power  spectrum  of  AT,  #(u,v),  is  then  evaluated  using  a 
2-dimensional  Fourier  transform.  Let  y  (and  y')  be  the  Fourier  transform  of 
y  (and  y'),  then 

y(u,v)  =  j  j  y(x,y)e-‘(“>«+''y)dxdy  (2.11) 

and  the  power  spectrun  associated  with  y  and  y' 

♦  =  y(u,v)  y'(u,v)*  (2.12) 

where  y'(u,v)*  is  the  complex  conjugate  of  y(u,v). 

The  power  spectrum  of  AT  is  then  obtained  from  ^|,t^'(u,v)  in  a  similar 
way,  as  in  the  spherical  case,  by  applying  the  inverse  linear  functionals 
related  to  y  and  y'  on  the  spectrum.  Then  the  2-dimenBional  covariance 
function,  K(x,y),  is  obtained  using  the  inverse  Fourier  transform. 

K(x,y)  =  4^  1  1  ♦(u,v)e*(«>‘+v»)dudv  (2.13) 

and  the  isotropic  covariance  function,  Kfs),  by  averaging  over  the  azimuth 
an 

K(s)  =  ^  I  K(x,y)do(  (2.14) 

o 

Using  this  procedure  the  integrations  in  eq.  (2.1)  are  calculated  as  a 
convolution  first  and  then  an  average  over  the  azimuth. 

If  the  power  spectrum,  4(u,v),  is  isotropic  (or  had  become  isotropic  by 
averaging  over  the  azimuth  in  the  frequency  domain)  the  isotropic  covariance 
function  K(s),  is  obtained  from  eq.  (2.13)  where  the  Fourier  transform  is 
reduced  to  an  inverse  Hankel  transform.  Consequently  the  isotropic  power 
spectrum  is  obtained  from  the  isotropic  covariance  function  using  a  Hankel 
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transform  (and  not  a  l-dimensional  Fourier  transform).  That  is 

♦  vcj)  -  I  K(s)Jo(scj)s  ds  (2.15) 

O 

and 

«• 

K(s)  =  I  ♦(«)Jo(6js)oj  d<j  (2.16) 

o 

where  =  u*+v*  and  is  the  Bessel  function  of  order  zero. 

For  details  cf.  Nash  and  Jordan  (1978),  Forsberg  (1984a)  and  Schwarz 

(1985). 

These  formulas  (2.11-2.16)  are  given  in  an  infinite  plane  and  the  spectrum, 
given  in  eq.  (2.12),  is  continuous,  but  it  becomes  discrete  in  the  local  case  if 
periodicity  is  assumed.  Then  the  integration  in  eq.  (2.11)  becomes  finite  and 
the  integration  in  eq.  (2.13)  and  (2.16)  reduces  to  a  summation.  The  discrete 
values  of  the  spectrum  appear  for  wavelengths  (xj-X|)/j  and  (ya-yi)/h  in 
each  direction,  when  j  and  k  are  positive  integers.  On  a  sphere  this 

corresponds  to  harmonic  degree  2»r»j/(xa-x,)  and  2rr*k/(y a-yi )  respectively. 

With  a  discrete  data  distribution  the  Discrete  Fourier  transform  is  used.  When 
data  are  arranged  in  a  regular  grid,  the  integrations  in  eq.  (2.11)  are 
calculated  as  sums  and  the  power  spectrum  becomes  periodic.  The  highest 
frequency  which  may  be  estimated  depends  on  the  spacing  of  the  data,  Ax  and 
^y,  since  the  smallest  wavelength  is  equal  to  two  times  the  spacing,  which  on 
a  sphere  corresponds  to  harmonic  degree  n/tx  and  n/Ay. 

For  details  on  the  application  of  the  Fourier  transform  cf.  Bracewell 

(1983),  and  Mesko  (1984). 

Eq.  (2.6)  and  eq.  (2.16)  express  the  covariance  function  in  a  spherical  and 
a  plane  approximation  respectively.  In  a  local  area  these  approximations 
converge  to  each  other  and  there  exists  a  link  between  the  degree  variances 
and  the  power  spectrum  (Porsberg,  1984a) 
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=  (i+l/2)l/(27TB)^((i+l/2)/R)  (2.17) 

where  R  is  the  distance  from  the  center  of  the  earth  to  the  plane.  The  left 
hand  side  expresses  the  degree  variances  in  eq.  (2.6)  on  a  sphere  with  radius 
K. 

Normalized  potential  degree  variances  are  obtained  as  dimensionless 
quantities  through  a  division  by  (GM/R)*. 

The  local  covariance  function  can  be  determined  in  two  ways.  The  first 
method  is  to  evaluate  eq.  (2.1)  using  eq.  (2.3)  or  eq.  (2.5).  The  other  method 
is  to  evaluate  eq.  (2.1)  (given  in  planar  coordinates)  using  the  Discrete 
Fourier  transform  and  an  azimuth-average.  The  advantage  of  the  second 
method  is  that  the  amount  of  computation  is  much  smaller  than  in  method  one, 
and  that  the  power  spectrum  is  obtained  during  the  computations.  The 
disadvantage  is  that  the  data  have  to  be  arranged  in  a  regular  grid. 

2.2  Estimation  of  a  local  empirical  covariance  function 

Before  a  set  of  observations  is  used  in  studies  of  the  gravity  field,  it  is 
important  that  non-gravimetric  signals  like  orbit  errors  and  sea  surface 
topography  in  altimeter  data,  are  removed  from  the  observations. 
Furthermore,  the  observations  must  be  associated  with  the  same  geodetic 
reference  system. 

Then  the  longwavelength  part  of  the  gravity  field  has  to  be  removed  in 
order  to  estimate  the  local  empirical  covariance  function.  If  this  is  done  using 
a  spherical  harmonic  approximation,  residual  observations,  ft  - 
L|(T)-L|(TN)'t’n|,  where  n,  is  the  noise  associated  with  the  i'th  observation, 
are  obtained.  Such  quantities  are  normally  used  in  studies  of  the  gravity 
field,  but  an  estimation  of  a  local  covariance  function  will  result  in  a 
covariance  function,  cov(y,y'),  that  is  affected  by  the  noise  of  the 
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observations  ,  and  errors  in  the  spherical  harmonic  approximation  (Kapp, 
1985).  That  is 

cov(y,y')  =  cov(y,y' )  +  cov(L(Tre#-TM) , L' {Tre^-T^) )  +  cov(n,n' ) 
where  cov<y,y'}  is  the  local  empirical  covariance  function,  C{ii), 
cov(L(Tref-Tn),  L' (Tp,f-T(|) )  is  the  covariance  function  associated  with  errors 
in  the  spherical  harmonic  approximation,  and  cov(n,n')  is  the  covariance 
function  associated  with  the  noise  of  the  observations.  The  remaining  terms 
are  assumed  to  be  zero.  Therefore,  the  estimated  covariance  funt;tion  is  nut 
the  local  empirical  covariance  function,  because  cov(L(Tr,f-T|*))  is  associated 
with  wavelengths  longer  than  the  extent  of  the  local  area.  Cov(L(Troj -Tk), 
L  (TrefTn))  is  indeed  a  part  of  the  covariance  function  that  should  be  used 
in  e.g.  collocation,  but  it  is  not  a  part  of  the  local  covariance  function  and 
can  never  be  estimated  from  observations  in  the  local  area.  Remaining  long 
wavelength  parts  of  the  gravity  field  will  interfere  with  the  result  and 
furthermore  cause  spectral  leakage  if  periodicity  if  assumed.  Consequently 
remaining  long  wavelength  parts  of  the  gravity  field  have  to  be  removed 
completely  from  the  observations.  Cov(n,n')  is  zero  for  ^  0,  if  the  noise  is 
assumed  to  be  uncorrelated.  For  cov(n,n')  is  equal  to  the  variance  of  the 
noise. 

In  a  local  area  remaining  long  wavelength  parts  of  the  gravity  field  may 
appear  as  a  bias.  In  Knudsen  (1987)  such  a  bias  was  removed  in  order  to 
center  the  observations  by  a  transformation  into  locally  best  fitting  referen«:e 
system.  However,  all  wavelengths  in  the  spherical  harmonic  approximation  may 
contain  errors.  Therefore  it  is  not  sufficient  to  remove  a  bias.  The 
procedure  must  be  able  to  consistently  remove  all  wavelengths  longer  than  the 
extent  of  the  local  area  from  several  kinds  of  gravity  field  related  quantities. 
A  method  that  fulfills  these  requirements  in  least  squares  collocation  using  a 
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covariance  function  that  is  designed  for  the  task.  Then  X|  =  L^iTr-pf-T^)  is 
estimated  from  7  =  {Lj(Tr,f-TH)+nj)  by 

xj  =  Ci  (C+D)~*  y  (2.19) 

and 

e*  =  -  Cj,(C+D)-*C«,  (2.20) 

where  C  is  a  matrix  containing  the  covariance  between  the  elements  in  y,  D  is 
a  matrix  containing  the  associated  error  covariances.  is  a  vector  holding 

covariance  values  between  the  estimate  and  y  and  is  the  variance  of 

the  quantities  x,  e  is  the  error  of  X|. 

In  a  spherical  harmonic  expansion  the  remaining  long  wavelength  part  of 
the  gravity  field  is  expressed  by  a  series  up  to  degree  and  order  N.  The 
covariance  function  between  quantities  related  to  this  part  of  the  gravity  field 
is  expressed  by  eq.  (2.6)  but  truncated  at  degree  N.  If  a  spherical  harmonic 
approximation  up  to  degree  and  order  M  has  been  subtracted,  where  M  is  less 
than  N,  the  degree  variances  from  the  global  empirical  covariance  function 
(see  Tscherning  and  Rapp  (1974))  are  used  from  degree  M-fl  to  N.  The 
quantities  LjiTraf-T),)  are  estimated  from  the  observations,  yj,  by  a 
convolution  with  a  filter  that  truncates  the  spectrum  at  harmonic  degree  N. 
This  has  to  be  done  in  the  space  domain,  because  the  long  wavelengths  are 
not  resolved  in  the  frequency  domain.  As  discussed  in  Jekeli  (1981),  a 
rectangular  filter  is  not  a  good  choice,  since  the  spectrum  of  this  filter  has 
relatively  large  sidelobes.  However,  no  expression  for  the  ideal  filter  is  given. 
In  the  plane  the  ideal  isotropic  filter  is  expressed  by  a  first  order 
Bessel  function.  A  reasonable  alternative  to  the  rectangular  filter  is  a  rec- 
tangular  sine  filter:  w(x,y)=sinc(-  N)  sinc(*'  N),  where  sinc(a)=sin(Tra)/(»ta) . 

Tf  7f 

Values  associated  with  several  kinds  of  gravity  field  related  quantities 
may  be  used  in  the  computations  simultaneously.  This  means  that  quantities. 


10 


Xj,  associated  with  different  kinds  of  gravity  field  related  quantities  are 
estimated  consistently.  It  is  also  possible  to  compute  corrections  to  the 
spherical  harmonic  expansion  (Tscherning,  1974). 

The  Fortran  program  GEOCOI,  (Tscherning,  1985)  can  be  used  for  the 
estimation  of  the  long  wavelength  part  of  the  gravity  field,  but  a  few  changes 
are  needed.  The  subroutine  COVAX  shall  be  called  with  LSUM  =  .TRUE.,  N2  is 
equal  to  the  degree  N,  and  HMAX  =  -1.0  (see  input  7  in  the  program). 
Furthermore  the  constant  A  (in  eq.(2.7))  shall  enter  CCI(8)  directly  without 
being  calculated  from  the  variance  of  gravity  anomalies.  In  the  subroutine 
the  array  SM  must  have  a  dimension  of  at  least  N. 

The  local  empirical  covariance  function  may  then  be  estimated,  when  the 
remaing  long  wavelength  parts  of  the  gravity  field  are  removed.  In  principle 
this  covariance  function  is  expressed  by  an  infinite  sum  of  a  series  of 
Legendre  polynomials,  but  in  practice  it  is  estimated  from  discrete 
observations  using  numerical  integration.  A  reasonable  result  can  only  be 
obtained  if  the  spectrum  tends  to  zero  and  the  spacing  between  the 
observations  is  so  small  that  aliasing  effects  caused  by  the  higher  frequencies 
are  negligible.  Since  the  gravity  field  is  unknown,  it  is  impossible  to 
determine  how  dense  the  observations  are  needed.  Therefore  the  spacing 
between  the  observations  has  to  be  determined  from  experiences  from  other 
investigations. 

A  global  empirical  covariance  function  was  estimated  by  Tscherning  and 
Rapp  (1974).  Then,  the  degree-variances  were  modelled  using  eq.  (2.7),  and 
the  factor  A  and  the  radius  of  the  Bjerhammar  sphere  were  adjusted  in  order 
to  fit  the  covariance  function  to  the  empirical  covariance  function.  The 
procedure  resulted  in  a  depth  to  the  Bjerhammar  sphere,  R-Ro,  of  1.22  km. 
Results  from  later  investigations  indicate  that  this  depth  should  be  larger. 
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For  the  Canadian  covariance  function  a  depth  of  3.35  km  was  found  by 
Schwarz  and  Lachapelle  (1980).  In  the  Faeroe  Islands  region  a  depth  of  3.17 
km  was  found  by  Knudsen  (1987).  In  Germany,  however,  a  depth  of  1.00  km 
was  found  by  Denker  and  Wenzel  (1987). 

Other  studies  have  been  carried  out,  where  the  decay  of  the  potential 
degree  variances  has  been  estimated  from  the  power  spectrum.  Such  studies 
indicate  (see  e.g.  Rapp,  1985,  and  Forsberg,  1987)  that  values  around  -3.6  are 
typical,  but  values  ranging  from  -3.3  to  -4.8  were  found  in  local  studies  of 
the  gravity  field  in  the  nordic  countries. 

In  order  to  compare  the  results  obtained  by  adjusting  Rg  and  the  decay 
of  the  potential  degree-variance,  the  decay  of  the  degree  variance  model 
multiplied  by  (Rg/R)**'*'^  was  computed.  Since  (Rg/R)**"*"*  tends  faster  to  zero 
than  a  polynomium,  the  decay  was  computed  in  interval  as  the  difference 
between  the  logarithms  of  the  degrees  variances  multiplied  by  (Rg/R)*''*’*  at 
harmonic  degree  180  and  1800.  With  a  depth  to  the  Bjerhammar  sphere  of  I.O, 
2.0,  3.0,  and  4.0  km.  a  decay  of  -3.2,  -3.4,  -3.6,  and  -3.8  was  found. 
Therefore,  a  depth  to  the  Bjerhammar  sphere  of  3.0  km  agrees  quite  well  with 
values  of  the  decay  around  -3.6.  However,  depths  around  1.5  km 
(correspKinding  to  decays  around  -3.3)  may  occur  as  well  as  larger  depths  to 
the  Bjerhammar  sphere. 

The  covariance  function  models  described  above  may  then  be  used  to 
evaluate  how  dense  the  observations  need  to  be  distributed  in  order  to 
estimate  the  local  empirical  covariance  function.  Degree-variances  associated 
with  different  kinds  of  obsevations  may  be  calculated  and  information  about 
how  fast  they  tend  to  zero  is  obtained.  Also  covariance  values,  Cf(,  may  be 
computed  from  different  harmonic  degrees  to  infinity  (Tscherning,  1976)  and 
information  about  the  magnitude  of  the  high  frequency  part  of  the  gravity 
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field  is  obtained.  Such  values  have  been  computed  for  jfeoid  heif^hts  and 
gravity  anomalies  with  a  depth  to  the  Bjerhammar  sphere  of  1.5  km  and  3.0 
kin.  In  this  case  the  gravity  field  above  harmonic  degree  180  is  assumed  to 
be  studied,  so  the  values  are  given  relative  to  the  values  at  degree  180  in  dR 
(101ogio(tfN/<^iao)  and  101ogio(0l|^/Ciao))  degree  360,  720,  1440,  and  2880 
(see  Table  1). 


(H-Ra) 

N 

_ _ 

_ Cm/Cibo _ 

(C.C) 

(Ag.dg) 

(C.O 

( Ag.Ag) 

1.5 

360 

-  9.17 

-  3.13 

~  6.45 

-  1.54 

720 

-18.82 

-  6.74 

-13.59 

-  3.70 

1440 

-29.26 

-11.16 

-21.73 

-  6.88 

2880 

-41.21 

-17.08 

-31.72 

-11.83 

3.0 

360 

-  9.54 

-  3.50 

-  7.03 

-  2.08 

720 

-19.92 

-  7.85 

-15.13 

-  5.27 

1440 

-31.84 

-13.74 

-25.09 

-10.65 

2880 

-46.73 

-22.60 

-38.43 

-18.36 

Table  1.  Values  of  degree-variances,  and  variances,  Cn,  in  dB 
relative  to  values  at  degree  180.  Values  associated  with 
geoid  heights,  and  gravity  anomalies,  (Ag.Ag),  are 

calculated  using  the  Tscherning/Rapp  model  4  and  a  depth  to 
the  BJerghammar  sphere,  R-Rg,  of  1.5  km  and  3.0  km.  The 
variances  are  calculated  as  the  sum  of  the  respective  degree- 
variances  from  harmonic  degree  N  to  infinity. 

First  of  all,  the  results  illustrate  that  the  degree-variances  associated 
with  geoid  heights,  tend  much  faster  to  zero  than  the  degree-varianciis 
associated  with  gravity  anomalies.  Furthermore,  the  degree-variances 
associated  with  a  depth  to  the  Bjerhammar  sphere  of  3.0  km  tend  to  go  faster 
to  zero  than  the  degree-variances  associated  with  a  depth  of  1.5  km.  The 
part  of  the  variance  that  is  located  above  the  different  harmonic  degrees  tend 
to  zero  in  a  manner  similar  to  the  degree-variances,  but  not  so  fast.  This 
shows  that  if  a  harmonic  degree  is  determined  so  the  degree-variances  have 
decreased  to  a  certain  level,  then  the  variance  above  this  harmonic  degree  has 
not  decreased  to  this  level.  Therefore  the  spacing  that  is  necessary  to 
resolve  the  spectrum  of  the  gravity  field  with  negligible  aliasing  effects 
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should  be  determined  considering  the  variance  values.  If  the  level  is  chosen 
to  be  between  -15  dB  and  -20  dB  the  variance  of  that  part  of  the  spectrum 
that  causes  the  aliasing  has  decreased  to  1%-3X.  Theni  values  associated  with 
a  depth  to  the  Bjerhaminar  sphere  of  3.0  km  for  a  harmonic  degree  of  720  is 
found  for  geoid  heights  and  2880  is  found  for  gravity  anomalies.  This  means 
that  geoid  heights  and  gravity  anomalies  are  needed  with  a  spacing  of  1/4* 
and  1/16*  respectively,  in  order  to  resolve  the  two  kinds  of  gravity  fields.  If 
a  depth  of  the  Bjerhammar  sphere  is  assumed  to  be  1.5  km  a  more  dense 
distribution  is  needed  for  both  kinds  of  observations. 

In  order  to  compute  covariance  values  using  eq.  (2.5),  the  local  area  is 
subdivided  into  small  cells  and  one  observation  in  each  cell  is  selected  in 
order  to  obtain  a  more  homogenously  distributed  data  set  with  a  spacing  close 
to  the  one  that  is  required.  Also,  data  outside  the  area  should  be  used. 
Periodicity  may  be  assumed,  but  then  problems  due  to  wavelengths  that  are 
not  periodic  may  arise.  As  discussed  in  Knudsen  (1987),  estimated  values 
should  not  be  filled  in  empty  cells. 

If  the  covariance  function  is  computed  using  the  discrete  Fourier 
transform  and  an  azimuth  average,  the  observations  have  to  be  gridded.  This 
can  be  done  using  least  square  collocation.  However,  faster  methods  like 
weighted  means  and  collocation  using  only  the  closed  observations  may  be 
used.  The  gridding  procedure  may  have  a  considerable  smoothing  effect 
(Knudsen,  1987),  which  can  be  diminished  by  using  a  much  more  densely 
distributed  set  of  observations  than  the  resulting  grid. 

Before  the  Fourier  transform  is  calculated,  it  is  advisable  to  apply  a 
cosine  tapered  window  in  order  to  avoid  spectral  leakage  caused  by 
wavelengths  that  are  not  periodic  in  the  local  area. 

A  cosine  taper  being  effective  on  K  points  from  each  border  is 
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A  cosine  taper  being  effective  on  K  points  from  each  border  is 


-  -  —  cosCn  ^7^)  for  k=l,...,K 
2  2  K 

1  for  k=(k+l N/2 


(2.21) 


where  N  is  the  number  of  points  in  each  direction. 

The  loss  of  power  may  be  re-established  through  a  multiplication  factor, 
so  agreement  with  the  variances  before  and  after  the  cosine  tapering  is 
obtained.  That  is 


yfj  =  c  y,j  IV,  wj  (2.22) 

where  c  is  the  nultiplication  factor 

c*  =  I  yfj  /  I  (y,j  w,  Wj)*  (2.23) 

tj  ‘J 

where  vT.j  “i"®  so-called  windowed  observations. 

Figure  1  and  2  show  the  effects  of  using  a  cosine  taper.  In  both  figures, 
the  power  spectra  (in  dB)  of  a  one-dimensional  discrete  Fourier  transform  of  a 
sequence  of  64  points  are  shown,  when  cosine  tapers  with  K=0  (do  nothing), 
K=8,  K=16,  or  K=32  are  used.  Figure  1  shows  the  results  where  a  cosine  with 
the  frequency  fo=16  was  used  as  sequence.  Then  the  Fourier  transforms  of 
the  tapers  centered  at  f=16  are  obtained,  since  the  multiplication  in  the  space 
domain  corresponds  to  a  convolution  in  the  frequency  domain.  These 
frequency  domain  impulse  responses  show  how  the  impulses  are  smoothened 
out,  when  K  increases.  The  shoulders,  however,  appearing  with  values  smaller 
than  -20  dB,  become  more  narrow.  With  K=32,  a  cosine  taper  similar  to  the 
Hanning  taper  is  obtained,  which  in  the  frequency  domain  has  the  values  1/4, 
1/2,  1/4  at  the  frequencies  -1,0,1  (see  eg.  Mesko,  1984). 
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Figure  1.  Effect  of  cosine  taper  with  K=0  (A),  K=8  (B),  K=16  (C), 
and  K=32  (D)  on  the  discrete  power  spectrum  of  a  cosine  function 
having  a  frequency  of  16.0  sampled  in  64  pqints. 

Figure  2  shows  how  spectral  leaka;(e  occurs  and  decreases  when  a  rosine 
taper  is  applies.  A  function  that  consists  of  10  cosines  with  the  frequencies: 
fg  =  15.6,... ,16.5  having  the  same  amplitude  (0.1)  was  used  as  sequence.  The 
power  of  this  sequence  is  expected  to  be  located  at  f=16  with  some  influence 
on  the  neighboring  frequencies.  If  no  cosine  taper  is  applied,  it  is  seen 
(figure  2A)  how  those  frequencies  are  resolved  in  the  discrete  spectrum.  Most 
of  the  power  (93.0%)  is  located  at  f=l5,  16,  and  17,  but  the  remaining  power 
(7.0%)  appears  at  all  frequencies  as  a  result  of  the  leakage  effects.  When  the 
cosine  tapers  are  applied  the  leakage  effects  decrease  and  the  parts  of  the 
power  that  are  located  as  f=l5,  16,  and  17  increase  to  93.8%,  96.0%,  and  99.4%. 
This  means  that  the  cosine  taper  with  K=32  (=N/2)  results  in  a  spectrum  with 
a  minimum  frequency  dispersion. 
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Figure  2.  Effect  of  cosine  taper  with  K=0  (A),  K=8  (B),  K=16  (C>, 
and  K=32  <D)  on  the  discrete  power  spectrum  of  a  sum  of  10  cosine 
functions  having  frequencies  of  15.6-16.5  sampled  in  64  points. 

The  estimaled  covariance  function  contains  both  the  local  empirical 
covariance  function  and  a  covariance  function  associated  with  the  noise.  In 
the  power  spectrum,  uncorrelated  noise  will  appear  as  a  constant  level  (white 
noise)  and  dominate  the  higher  frequencies,  where  the  power  of  the  gravity 
field  is  small.  The  noise  in  itself  is  unknown,  but  the  variance  of  the  noise 
may  be  estimated  if  reliable  noise  terms  are  associated  with  the  observations. 
Another  possibility  is  to  estimate  the  noise  level  in  the  power  spectrum.  Then 
the  estimated  variance  of  the  noise  is  subtracted  from  the  estimated 
covariance  function  and  the  local  empirical  covariance  functions  associated 
with  the  earth’s  gravity  field  is  obtained. 
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2.3  The  adjustment  of  a  covariance  function  model 

The  adjustment  of  a  covariance  function  model  is  an  important  step  when 
the  local  empirircal  covariance  function  is  used  in  estimations  of  gravity  field 
related  quantities.  In  Knudsen  (1987)  the  Tscherning/Rapp  model  was  fitted  to 
the  empirical  covariance  values  by  adjusting  the  depth  to  the  Bjerhammar 
sphere,  (R-Rg),  the  factor  A,  and  a  scale  factor,  a,  associated  with  the  error 
degree  variances  using  an  iterative  least  squares  inversion  procedure.  In 
each  iteration,  the  adjustment  of  the  parameters  Xg  were  calculated  as 

X  -  x„  =  (ATC-*A  +  ATCy-‘  (y-y„)  (2.24) 

where 

X  is  the  adjusted  paraaeter 

y  is  the  eapirical  covariance  values 

Vo  is  the  values  from  the  nodel  using  x^ 

A  is  the  Jacobian  natrix  {dy|/9xj} 

Cy  is  the  error  covariance  matrix  of  y 
C*  is  the  covariance  matrix  of  (x-Xj,) 

The  ablility  of  the  model  to  describe  the  empirical  values,  or  the  fitness,  was 
measured  by  the  dimensionless  Q  value,  where 

0*  =  (y-Yo)^  Cy‘(y-yo)  (2.25) 

where  n  is  the  number  of  data  and  m  is  the  number  of  parameters  (=3). 

In  practice  the  adjustments  were  calculated  relatively  as  dimensionless 
quantities  by  multiplying  each  column  in  matrix  A  by  the  associated  parameter. 
Consequently,  both  the  rows  and  columns  in  matrix  were  divided  by  the 
associated  parameters  obtaining  relative  apriori  variances. 

The  error  covariance  matrix,  Cj,  was  assumed  to  be  diagonal  and  contained 
the  square  of  some  empirical  error  estimates.  These  errors  were  calculated  in 
order  to  evaluate  the  accuracy  of  numerical  integration  in  eq.  (2.5)  and 
depends  on  the  variances,  Cg  and  Cg  of  the  observations  y  and  y',  the  size 
of  the  local  area,  the  size  of  the  cells.  At  and  AX,  and  the  actual  number  of 
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products  that  would  appear  if  one  observation  is  located  in  each  and  every 
cell.  The  error,  errj,,  associated  with  the  k'th  covarianca  value  is  then 


err^ 


where 


3x 

Nk 


C--  y~c„c: 


n^- 


)  (X.-X, ) 

2  A*AX 


for  k=0 
for  k  ^  0 


12.2G1 


The  method  was  tested  in  the  Faeroe  Islands  region  using  throe 
combinations  of  empirical  covariance  values  as  input.  The  results  showed  that 
the  scale  factor,  a,  was  not  well  determined  from  gravity  anomaly  covariance 
values,  and  the  depth  to  the  Bjerhammar  sphere  was  not  well  determined  from 
geoid  height  covariance  values  only.  A  combination  of  these  two  kinds  of 
covariance  values  resulted  in  a  well  determined  model:  (R-Rb)=3.17  *0.34  km, 
A=889  *47  10’m*/s*,  and  a=0.21  *0.04.  Changes  in  the  initial  model  were  not 
found  to  have  influence  on  the  results. 

If  agreement  between  the  extent  of  the  area  and  the  degree,  which  the 
gravity  field  have  been  removed  up  to,  exists,  the  scale  factor  should  be  zero. 
This  was  not  the  case  in  the  Faeroe  Islands  region.  The  scale  factor  was  not 
zero  because  the  remaining  longwavelength  parts  were  not  completely  removed 
by  the  centering  that  was  used.  Furthermore,  the  area  was  larger  than  the 
wavelengths  that  were  removed.  Suppose  that  a  spherical  harmonic  expansion 
up  to  degree  M  had  been  subtracted  from  the  observations,  and  that 
remaining  long  wavelength  parts  have  been  removed  up  to  degree  N  (N^M) 
using  the  procedure  described  in  section  2.2.  N  corresponds  to  the  size  of 
area.  Then  the  degree  variances  up  to  degree  N  are  zero.  The  degree 
variances  between  degree  N  and  M  are  a  part  of  the  local  covariance  function 
and  they  may  be  modelled  by  the  error  degree  variances  that  are  associated 
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and  they  may  be  modelled  by  the  error  decree  variances  that  are  associated 
with  the  spherical  harmonic  expansion,  multiplied  by  a  scale  factor. 

The  use  of  a  local  covariance  function  model  from  degree  N-t-1  and  a  global 
covariance  function  model  up  to  degree  N,  when  quantities  like  9  is  used,  is 
discussed  in  section  2.2.  The  local  covariance  function  is  also  of  interest  in 
estimation  methods  like  point  mass  modelling.  Porsberg  (1984a)  describes  the 
relations  between  least  squares  collocation  and  point  mass  modelling  and  how 
the  depth  to  the  point  masses  affects  the  shape  of  the  covariance  function.  A 
covariance  function  with  an  asymptotic  decay  of  the  potential  degree  variances 
of  -3  on  the  surface  of  a  Bjerhammar  sphere  (like  the  Tscherning/Rapp  model) 
is  implicitly  obtained  when  a  certain  distribution  of  point  masses  is  located  at 
a  depth,  which  is  two  times  the  depth  to  the  Bjerhammar  sphere. 

3.  The  Use  of  Satellite  Altimeter  Data  for  Estimation  of  Local  Empirical 

Covariance  Functions 

In  this  chapter  three  local  empirical  covariance  functions  are  estimated 
from  locally  crossover  adjusted  Seasat/Cieos>3  altimeter  data.  The  purpose  of 
a  local  crossover  adjustment  is  described  in  a  brief  discussion  on  the  use  of 
satellite  altimetry  and  the  results  from  the  adjustments  of  the  altimeter  data 
are  evaluated.  Covariance  functions  associated  with  the  gravity  field  above 
harmonic  degree  180  are  estimated  and  the  spherical  harmonic  expansion  OSU81 
(Rapp,  1981)  are  used  as  reference.  The  effects  of  remaining  long  wavelength 
parts  on  the  estimation  of  empirical  covariance  functions  are  studied  by 
comparing  covariance  functions  estimated  before  and  after  those  remaining 
long  wavelength  parts  of  the  gravity  field  were  removed  from  the 
observations. 

In  order  to  study  the  variability  of  the  gravity  field  the  three  local 
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areas  were  selected  with  different  characteristics.  Area  1  (The  New  Enjtland 
Seamount  Area:  38*<4<40*,  295*<X<297*)  was  selected  as  a  ’mild'  area,  Area  2 
(34*<4<36*,  294*<X<296*)  was  selected  as  a  ’smooth'  area,  and  Area  3  (The 
Bermuda  Area:  3l*<4<33*,  294*<X<296*)  was  selected  as  a  ’rou^h’  area. 

3.1  The  Crossover  Adjustment 

Prom  satellite  altimeter  data  sea  surface  heii^hts  are  derived  usin^  the 
ellipsoidal  heights  of  the  satellite.  These  sea  surface  heights  may  have  been 
corrected  for  a  number  effects  (instrumental  and  atmospherical  effects  and  sea 
state  related  bias).  Prom  such  instantaneous  sea  surface  heights  mean  sea 
surface  heights  are  obtained  by  subtracting  variations  in  the  sea  surface 
heights  due  to  tide  and  variations  in  the  atmospheric  pressure.  After  a 
removal  of  the  sea  surface  topography  observations  of  geoid  heights  are 
obtained.  The  accuracy  of  such  observations  depends  on  the  quality  of  the 
models  that  have  been  used  in  the  corrections  of  the  altimeter  data. 
Furthermore  the  observations  contain  unmodelled  phenomena  like  effects  from 
rain,  clouds,  and  changes  in  ocean  currents.  A  geoid  height  observation 
derived  from  an  altimetric  observation  may  therefore  be  described  by 

h  =  f  +  Ahc  +  Aht  +  n  (3.1) 

where  C  is  the  geoid  height  and  the  non-gravimetric  signal  is  divided  into  a 
constant  part,  4h(.,  and  a  time  varying  part,  Ah^.  n  is  the  noise  of  the 
observation. 

Prom  a  set  of  observations,  (ht),  located  in  a  local  area,  a  quantity  x  may 
be  estimated  by  least  squares  collocation 

X  =  C^h  (Chh  +  D)-‘  {h,}  (3.2) 

where  C^i,  is  a  vector  containing  covariance  values  between  x  and  the 
observations  and  is  a  matrix  containing  covariance  values  between  the 
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observations.  D  is  the  covariance  matrix  associated  with  the  noise  of  the 
observations.  The  covariance  values  associated  with  the  observations  depend 
on  the  covariances  between  the  the  Ahc  and  the  Ah^  terms  in  eq.  (3.1),  and 
the  quantity  x  may  be  associated  with  at  least  one  of  those  terms.  A 
separation  of  the  geoid  and  the  constant  part  of  the  non-gravimetric  signal 
will  only  be  successful!  if  the  spectral  characteristics  of  the  two  signals  are 
differentt  and  a  separation  of  the  time  varying  components  will  only  be 
successful!  if  observations,  where  the  constant  parts  of  the  signal  (including 
the  geoid)  are  strongly  correlated,  are  available.  The  first  criterion  may  be 
fulfilled  if  the  constant  part  of  the  non-gravimetric  signal  consists  of  long 
wavelengths  and  those  wavelengths  have  been  subtracted  from  the  gravity 
field  related  part  of  the  signal.  The  second  criterion  may  be  fulfilled  if 
repeat  observations  are  used. 

The  need  of  repeat  observations,  when  quantities  are  estimated  using  eq. 
(3.2)  from  observations  containing  time  varying  components,  may  result  in 
very  large  equation  systems  and  problems  in  solving  them.  It  is  therefore 
desirable  to  remove  the  time  varying  components  from  the  observations  so 
repeat  observations  no  longer  are  needed.  The  problem  is  then  to  estimate 
the  time  varying  components,  because  the  repeat  observations  still  are  needed 
in  this  step.  In  this  case  a  solution  is  obtained  by  using  the  differences 
between  the  repeat  observations.  That  is  using  eq.  (3.2)  and  pairs  of 
observations,  h|  and  h'|,  where  h^  and  h*|  are  located  at  the  same  point  on 
either  colinear  arcs  or  crossing  arcs 

f(t)  =  C|  at  aT~‘  (Chh+U)-‘  A-‘  A  {h,,  h',} 

=  C{  AT  (A(Chh+D)AT)-‘  {d,} 

=  Cl  AT  (C’+D’)-*  {d,}  (3.3) 
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where  f(t)  is  the  estimate  of  Ah^  in  eq.(3.1).  {d()=A{h|,  h’|}  is  a  set  of 

differences  dt=h|-h’4.  If  the  time  varying  components  are  uncorreiated  with 
the  constant  terms  (f  and  Sh^)  contains  the  covariance  values  cov(4ht(t), 

Ah^t)  end  cov(Aht(t),  Ah'^i)  end  C  contains  the  covariance  values  cov(Ah^i, 
Ahtj)  +  cov(Ah’tif  Ah'tj)  -  cov(Ahtn  Ah’^j)  -  cov(Ah'tjr  Ah^j).  D’  contains  the 
error  variances  of  the  differences  d|.  The  error  of  the  estimated  value  is 
expressed  by 

err»(t)  =  c  ~  Cj  (C’  +  D’)“*  A  (3.4) 

where  c  is  the  variance  of  the  time  varying  components. 

After  the  time  varying  components  are  estimated  the  altimeter  data  are 
corrected  and  crossover  adjusted  observations  are  obtained.  That  is 
h,  =  h  -  f(t) 

=  (  +  Ahg  +  n,  (3.5) 

where  n,  is  the  noiae  of  the  crossover  adjusted  observation. 

The  method  and  how  the  covariance  values  may  be  computed  are  discussed 
in  Knudsen  (1987a).  Tests  in  The  Faeroe  Islands  Region  using  the  adjusted 
Seasat  altimeter  data  (Liang,  1983)  showed  that  remaining  time  varying 
components  were  successfully  removed.  It  was  assumed  that  orbit  related 
errors  had  been  removed  and  that  the  main  part  of  the  remaining  variations 
was  caused  by  an  inaccurate  ocean  tide  model  and  the  unmodelled  phenomena 
mentioned  above.  The  covariance  values  where  calculated  using  a  gaussian 
function  and  the  along  track  distance  between  the  crossover  points.  As 
correlation  distance  1000  km  was  used. 

The  potential  of  this  method  is  that  altimeter  data  in  local  areas  with 
inaccurate  ocean  tide  models  (coastal  areas  like  the  Mediterranean)  can  be 
adjusted  more  accurately  than  in  conventional  bias  or  bias/tilt  adjustments. 
On  the  other  hand  it  is  not  felt  that  the  method  is  suitable  for  orbit  error 
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estimationsi  since  orbit  errors  are  geographically  correlated  (see  eg.  Kngelis, 
1987).  Further  studies  of  time  varying  components  and  their  correlations  are 
needed  and  can  for  example  be  based  on  the  extensive  repeat  information  made 
available  through  the  Geosat  mission. 

3.2  Results  from  the  Crossover  Adjustment  of  Seasat/Geos-S  altimeter  data. 

For  the  estimation  of  the  local  empirical  covariance  functions  observations 
from  the  merged  Seasat/Geos-3  data  set  (Liang,  1983)  were  used.  All 
observations  in  each  2*  x  2*  area  plua  a  2*  border  zone  were  retrieved.  It 
resulted  in  three  data  sets  each  containing  35000  (approximately)  observations. 
Since  the  observations  were  distributed  far  more  densely  than  needed,  it  was 
decided  to  thin  out  the  data.  During  this  process  a  cell  size  of  1/6*  by  1/6* 
was  chosen  for  the  estimation  of  the  empirical  covariance  functions  (explained 
in  the  next  section).  Then  a  subset  of  arcs  were  selected,  so  the  observations 
associated  with  those  arcs  would  provide  at  least  one  observations  in  each 
cell.  This  was  done  in  each  area  and  resulted  in  data  sets  each  containing 
12000  observations  approximately. 

Crossover  discrepancies  were  computed  as  the  difference  in  height 
between  ascending  and  descending  tracks.  The  heights  and  the  position  were 
calculated  by  a  linear  interpolation  between  the  four  neighboring  observations. 
The  crossover  adjustments  were  carried  out  using  the  method  described  in 
section  3.1.  It  was  assumed  that  orbit  errors  had  been  removed  and  a 
covariance  function  similar  to  the  one  used  in  the  Faeroe  Islands  Region  was 
used.  Since  the  number  of  crossover  discrepancies  was  large,  the  estimation 
of  the  time  varying  components  was  based  on  crossover  discrepancies  between 
Seasat  arcs,  and  crossover  discrepancies  between  Seasat  and  Geos-3  arcs  only. 
RMS  values  of  the  discrepancies  were  computed  before  and  after  the 
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adjustments.  These  values  are  shown  in  Table  2. 

The  results  from  the  crossover  adjustments  (Table  2)  show  that 
discrepancies  between  Seasat  arcs  generally  are  much  smaller  than 
discrepancies  between  Geo8-3  arcs.  The  reason  for  that  is  that  the  Seasat 
altimeter  data  are  much  more  accurate  than  the  Geos-3  altimeter  data.  The 
results  show  furthermore  that  the  discrepancies  after  the  adjustments  are 
largest  in  Areas  1  and  smallest  in  Area  3.  It  was  also  expected,  since  Area  1 
is  located  in  the  center  of  an  area  with  large  sea  surface  variations  due  to 
changes  in  the  Gulf  Stream  (Menard,  1983).  Those  variations  may  contain  a 
rather  short  wavelength  signature,  which  is  not  completely  modelled  in  the 
adjustment. 


Area 

Type 

Number 

fiMS» 

RMS*> 

1 

S/S 

0.41  ■ 

0. 15  D 

S/G 

0.52  - 

0.34  - 

G/G 

1737 

0.60  - 

0.58  - 

2 

S/S 

106 

0.25  ■ 

0.13  n 

S/G 

964 

0.48  - 

0.31  - 

G/G 

2084 

0.60  - 

0.54  - 

3 

S/S 

111 

0. 13  n 

0.09  n 

S/G 

1064 

0.37  - 

0.27  - 

G/G 

2322 

0.51  - 

0.46  - 

Table  2.  Results  of  the  crossover  adjustments  in  Area  1,  2,  and  3  extended 
with  a  border  zone  of  2  degrees.  The  number  and  RMS  values  before  (a)  and 
after  (b)  the  adjustments  of  discrepancies  between  Seasat  arcs  (S/S),  between 
Seasat  and  Geos-3  arcs(S/G),  and  between  Geos-3  arcs  (G/G). 


Crossover  adjusted  altimeter  data  were  obtained  using  eq.  (3.5)  and  their 
noise  terms  were  evaluated.  The  purpose  of  an  evaluation  of  the  noise  terms 
is  to  determine  whether  they  are  appropriate  geoid  observation  noise  terms  or 
not.  As  mentioned  in  section  3.1  the  accuracy  depends  on  the  quality  of  the 
terms  associated  with  the  data  do  not  take  this  into  account.  Prom  studies  of 
repeat  tracks  by  Marks  and  Sailor  (1986)  typical  geoid  noise  terms  of  8  cm  for 
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Seasat  and  28  cm  for  Geos-3  altimeter  data  were  found.  These  values  are 
approximately  40%  higher  than  the  typical  errors  of  6  cm  and  20  cm  that  are 
assigned  Seasat  and  Geos-3  altimeter  data  respectively.  In  this  case  the  noise 
terms  are  evaluated  using  RMS  values  of  the  crossover  discrepancies  divided 
by  the  errors  of  the  discrepancies.  If  the  numbers  are  larger  than  one,  time 
varying  components  are  still  present.  Such  remaining  time  varying 
components  are  correlated  along  track,  but  across  track  they  may  be  treated 
as  uncorrelated  errors.  If  the  correlations  along  track  are  small,  the  total 
errors  may  be  assumed  to  be  uncorrelated.  Consequently  the  noise  terms  can 
be  modified,  so  they  agree  with  the  magnitude  of  the  crossover  discrepancies. 
This  was  done  in  each  area  by  computing  scale  factor  for  the  noise  terms  of 
Seasat  and  Geo8-3  respectively.  Then  the  noise  terms  were  multiplied  by  the 
respective  scale  factors  and  noise  terms  associated  with  geoid  height 
observations  were  obtained.  The  scale  factors  and  RMS  values  of  the  modified 
noise  terms  are  shown  in  Table  3. 


Area 

Seasat 

Geos-3 

factor 

RMS 

factor 

RMS 

1 

1.58 

0.11  m 

2.03 

0.41  m 

2 

1.46 

!.%• 

3 

1.46 

Table  3.  Scale  factors  and  RMS  values  of  the  modified  noise  terms  for  Seasat 
and  Geos-3  in  Area  1,  2,  and  3. 

3.3  The  Estimation  of  the  Local  Empirical  Covariance  Functions. 

It  was  decided  to  estimate  local  empirical  covariance  functions  associated 
with  harmonic  degrees  greater  than  180,  which  roughly  corresponds  to 
wavelengths  shorter  than  2  degrees,  and  use  the  space  domain  method  with 
observations  selected  in  cells  covering  the  areas.  The  size  of  the  local  areas 


26 


areas  should  be  taken  into  account  within  distances  of  2*.  Then  covariance 
values  are  calculated  usin^:'  eq.  (2.5)  by  forming  averages  of  products  between 
observations  in  the  local  2*  by  2*  areas  and  observations  in  the  6*  by  6* 
areas.  It  corresponds  to  a  typical  estimation  situation,  where  quantities  in  a 
local  area  are  estimated  from  observations  in  the  local  area  and  a  border  zone. 
The  variances  of  the  modified  noise  terms  are  also  calculated  and  subtracted 
from  the  variances  of  the  observations. 

The  cell  size  was  determined  from  the  results  in  section  2.2  and  from 
results  by  Marks  and  Sailor  (1986).  The  results  from  section  2.2  showed  that 
a  cell  size  of  1/4*  would  be  sufficient  to  resolve  geoid  heights.  Marks  and 
Sailor,  however,  found  that  the  short  wavelength  resolution  limit  is  about  32 
km  for  Seasat  and  about  60  km  for  Geos-3.  These  wavelengths  corresponds 
roughly  to  sample  spacings  of  1/7*  and  1/4*  (half  wavelength).  Since  a 
combination  of  Seasat  and  Geo8-3  data  are  used,  a  cell  size  of  1/6*  was 
chosen.  Then  each  area  was  subdivided  into  cells  of  1/6*  by  1/6*  and  one 
observation  in  each  cell  were  selected  from  the  locally  crossover  adjusted 
altimeter  data. 

The  removal  of  the  information  content  of  wavelengths  longer  than  the 
extent  of  the  local  areas  was  in  the  first  place  done  by  subtracting  the 
contribution  from  the  spherical  harmonic  expansion  OSU81  up  to  degree  and 
order  180.  Then  empirical  covariance  functions  were  estimated.  This  was 
done  in  order  to  study  the  effects  of  remaining  long  wavelength  parts. 

Then  remaining  long  wavelength  parts  were  estimated  using  the  method 
described  in  section  2.2.  Mean  values  were  computed  using  a  rectangular  sine 
filter  on  a  1/2*  by  1/2*  grid  covering  each  6*  by  6*  area.  From  these  mean 
values  remaining  long  wavelength  parts  in  the  observations  were  estimated  by 
least  squares  collocation  (eq.  (2.19))  and  a  covariance  function  truncated  at 
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leutit  squares  collocation  (eq.  (2.19))  and  a  covariance  function  truncated  at 
degree  180.  The  variances  of  the  mean  values  were  smaller  than  the  error 
variance  of  OSU81  (1.15  m^),  so  constant  non-gravimetric  effect  were  assumed 
to  be  small  and  included  in  the  estimated  values,  since  those  effects  mainly 
are  due  to  sea  surface  topography  of  wavelengths  corresponding  to  harmonic 
degrees  smaller  than  180.  After  a  removal  of  the  estimated  values  from  the 
crossover  adjusted  altimeter  data,  observations  of  geoid  heights  associated 
with  harmonic  degrees  greater  than  180  were  obtained.  Then  local  empirical 
covariance  functions  were  estimated  from  these  space  domain  highpass  filtered 
observations. 

The  estimated  geoid  height  variances  in  the  three  local  areas  were  0.256 
m*,  0.372  m*,  and  1.899  m*  respectively,  before  the  remaining  long  wavelength 
parts  (after  the  subtraction  of  the  0SU81  field)  were  removed.  (Noise 
variances  of  0.069  m*,  0.090  m*,  and  0.033  m*  respectively  have  been 

subtracted.)  Using  the  highpass  filtered  observations  these  numbers  reduced 
to  0.213  m’,  0.112  m^,  and  1.217  m^.  The  variances  associated  with  the 
remaining  long  wavelength  parts  are  then  0.043  m’,  0.260  m*,  and  0.682  m^  or 
20%,  232X,  and  56X  relative  to  the  variances  of  the  highpass  filtered 

observations.  The  effects  of  the  remaining  long  wavelength  parts  on  the 
estimation  of  the  geoid  height  variances  are  therefore  considerable.  Especially 
the  results  obtained  in  Area  2,  where  the  magnitude  of  the  remaining  long 
wavelength  parts  is  more  than  twice  as  big  as  the  magnitude  of  the  highpass 
filtered  observations,  show  the  importance  of  removing  those  long  wavelength 
parts. 

The  local  empirical  covariance  functions  estimated  in  Area  1,  Area  2,  and 
Area  3  using  the  selected  observations  before  and  after  the  removal  of  the 
remaining  long  wavelength  parts  are  shown  in  Figure  3-5.  Furthermore  power 
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spectra  calculated  from  1/6*  by  1/6*  mean  values  using  a  cosine  taper  (eq. 
(2.21  >)  with  K=6  are  shown.  The  purpose  of  calculating  power  spectra  was  to 
evaluate  the  spectral  content  at  the  lower  frequencies,  so  no  attempts  to 
correct  for  effects  due  to  smoothing  and  noise  were  made.  The  local  areas 
were  extended  to  4*  by  4*  in  order  to  estimate  the  1/4  cycle/degret!  power 
also.  The  effects  from  the  remaining  long  wavelength  parts  are  also  seen  by 
comparing  the  two  covariance  functions  in  each  area.  Before  tht;  highpa.ss 
filtering  the  covariance  functions  in  Area  1  and  Area  2  have  their  first  zero 
crossing  at  lags  of  1.75  deg.  and  more  than  2.00  deg.  respectively.  After  the 
highpass  filtering  the  first  zero  crossings  occurs  at  lags  close  to  0.5  deg. 
The  power  spectra  show  the  different  spectral  contents.  Before  the  highpass 
filtering  the  observations  contain  large  signals  associated  with  frequencies  1/4 
and  1/2  cycles/degree,  which  efficiently  are  removed  from  the  highpass 
filtered  observations.  The  information  contents  of  wavelengths  longer  than 
the  extents  of  the  local  areas  have  therefore  been  removed  from  the  highpass 
filtered  observations.  Consequently  the  covariance  functions  calculated  from 
those  highpass  filtered  observations  are  estimates  of  the  local  empirical 
covariance  functions  associated  with  the  gravity  field  above  harmonic  degree 
180. 
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Figure  3.  Upper:  Local  empirical  covariance  functions  calculated 
in  Area  1  from  altimeter  data  minus  OSU81  (o)  and  from  high-pass 
filtered  altimeter  data  (+)  (units  in  m* ) .  Lower:  Power  spectra 
calculated  in  Area  1  from  lO'xlO'  mean  values  of  altimeter  data 
minus  OSU81  (o),  high-pass  filtered  altimeter  data  (+),  and  their 
differences  ( x ) . 
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Figure  4.  Upper:  Local  empirical  covariance  functions  calculated 
in  Area  2  from  altimeter  data  minus  OSU81  (o)  and  from  high-pass 
filtered  altimeter  data  (  +  )  (units  in  m*  ) .  Lower:  Po  .er  spectra 
calculated  in  Area  2  from  lO'xlO*  mean  values  of  altimeter  data 
inus  OSU81  (o),  high-pass  filtered  altimeter  data  (+),  and  their 
1 f f erences  ( x ) . 
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Figure  5.  Upper:  Local  empirical  covariance  functions  calculated 
xn  Area  3  from  altimeter  data  minus  0SU81  (o)  and  from  high-pass 
filtered  altimeter  data  (+)  (units  in  ro*  ) .  Lower:  Power  spectra 
calculated  in  Area  3  from  lO'xlO'  mean  values  of  altimeter  data 
minus  0SU81  (o),  high-pass  filtered  altimeter  data  (-*-),  and  their 
differences  ( x  )  . 
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t.  Local  Empirical  Covariance  Functions  from  Residual  Terrain  Reduced 

Altimeter  Data. 

The  effect  from  the  topography  is  a  considerable  part  of  the  gravity  field 
and,  as  far  as  the  topography  is  known,  this  effect  can  be  used  when  the 
gravity  field  is  modelled.  This  may  be  done  by  computing  gravity  field 
related  quantities  from  the  topography  and  subtract  the  terrain  effects  from 
the  observations.  Then  terrain  reduced  observations  are  used  in  the 
estimations  and  the  results  are  obtained  by  adding  the  estimated  quantities 
and  the  terrain  effects.  If  a  spherical  harmonic  expansion  is  used  as 
reference,  terrain  effects  associated  with  the  same  wavelengths  as  the 
harmonic  expansion  are  removed  implicitly.  Then  residual  terrain  effects  are 
used.  Methods  for  computing  residual  terrain  effects  from  a  digital  terrain 
model  using  rectangular  prisms  or  Fourier  techniques  are  desribed  in  e.g. 
Forsberg  (1984)  and  Porsberg  (1985). 

The  effect  of  using  topographic  information  is  that  the  magnitude  of  the 
unknown  parts  of  the  gravity  field  becomes  smaller.  Furthermore  strongly 
varying  gravity  fields  in  mountainous  areas  reduce  to  more  smooth  gravity 
fields.  These  effects  are  described  in  Porsberg  (1986),  where  results  from  a 
study  of  the  spectral  properties  of  the  gravity  field  in  the  Nordic  countries 
are  evaluated.  RMS  values  of  gravity  anomaly  observations  relative  to  GPM-2 
(Wenzel,  1985)  to  degree  and  order  180  were  calculated  in  38  2*  x  4*  blocks 
before  and  after  residual  terrain  effects  were  removed  from  the  observations. 
In  Figure  6  the  RMS  values  of  the  observations  relative  to  GPM-2  (FA: 
Residual  Free-air  Anomalies)  are  plotted  against  RMS  values  of  the  residual 
terrain  reduced  observations  relative  to  GPM-2  (BA:  Residual  Bouguer 

corrected  Anomalies).  Furthermore  the  distribution  of  the  two  types  of  RMS 
values  are  shown.  The  RMS  values  of  the  residual  free-air  anomalies  range 
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from  8  mgal  to  59  mgal.  Typical  values  of  10-15  mgal  were  found  in  lowland 
areas  in  Denmark  and  Finland,  and  typical  values  of  30-45  mgal  were  found  in 
mountainous  areas  in  Norway.  The  RMS  values  of  the  residual  Bouguer 
corrected  observations  range  from  7  mgal  to  21  mgal.  The  mean  value  of  the 
RMS  values  dropped  from  19.3  mgal  to  12.3  mgal,  and  the  standard  deviation  of 
the  RMS  values  dropped  from  12.5  mgal  to  3.3  mgal,  when  the  residual  terrain 
effect  were  removed.  Areas  with  smooth  gravity  fields  and  smooth  topography 
did  not  change  much,  but  the  rough  areas  became  mxach  smoother.  The  most 
dramatic  change  was  found  in  a  Norweigian  Block,  where  the  RMS  value 
dropped  from  59  mgal  to  14  mgal.  The  gravity  field  in  the  Nordic  countries 
has  thereby  become  much  more  homogeneous,  and  the  magnitude  of  unknown 
parts  of  the  gravity  field  did  decrease. 

In  this  chapter  a  method  for  computing  residual  terrain  effects  from  an 
isostatic  earth  model  is  described.  The  method  is  used  in  the  three  local 
areas  described  in  the  previous  chapter  and  local  empirical  covariance 
functions  are  estimated  from  residual  terrain  reduced  altimeter  data. 
Furthermore  correlations  between  the  altimetric  and  bathymetric  geoid 
undulations  are  evaluated. 

4.1  Calculation  of  Residual  Geoid  Undulations  from  an  Isostatic  Earth  Model. 

Since  the  earth  is  believed  to  be  about  90X  isostically  compensated,  it  was 
decided  that  the  computation  of  geoid  undulations  from  the  topography  should 
take  the  isostasy  into  account.  Therefore  an  isostatic  earth  model  was  needed. 
A  simple  (and  highly  idealized)  model  is  the  Airy-Heiskanen  model  (Heiskanen 
and  Moritz,  1967).  This  model  is  based  on  a  floating  theory,  where  the 
topography  is  compensated  at  a  depth  of  30  km  by  root/antiroot  formations  at 
the  crust/mantle  interface.  Seismic  results  indicate  that  the  depth  to  the 
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crust/mantle  interface  is  correlated  with  the  topography,  which  may  justify 
the  principle  of  an  isostatic  compensation  at  the  crust/mantle  interface.  Thu 
Airy-Heiskanen  model,  however,  assumes  a  strictly  local  compensation,  which 
more  likely  has  a  regional  character  due  to  the  elasticity  of  earth  materials. 
The  Vening  Meinesz  model  assumes  such  a  regional  compensation  and  is  a 
"smeared  out"  version  of  the  Airy-Heiskanen  model.  Even  though  the  Vening 
Meinesz  model  is  more  realistic  than  the  Airy-Heiskanen  model,  the 
compensation  mechanisms  are  much  more  complicated.  Experiments  by  Lewis 
and  Dorman  (1970)  suggests  that  the  compensation  takes  place  at  different 
levels  depending  on  the  wavelengths  of  the  topography. 

It  was  decided  to  use  an  isostatic  Airy-Heiskanen  model  assuming  that  a 
gravity  field  similar  to  the  gravity  field  from  a  Vening  Meinesz  regional  model 
is  obtained  by  increasing  the  depth  of  compensation.  An  increase  of  the 
depth  of  compensation  corresponds  to  a  smoothing  of  the  crust/mantle 
topography  by  an  upward  continuation.  Then  the  depth  to  the  crust/mantle 
interface  is  given  by 

d'  =  Dc  -  o(d  (4.i) 

where  d  is  the  bathymetric  depth,  is  the  depth  of  compensation,  and 
(x=Ap^/Apc’  is  density  constant  associasted  with  the  topography  and  is 

the  crust/mantle  density  contrast.  As  densities  of  ocean  water,  crust,  and 
mantle  values  of  1.03  g/cm’,  2.67  g/cm’,  and  3.27  g/cm*  were  assumed.  Then 
A^i=1.64  g/cm’,  g/cm^,  and  o(=Z.73. 

The  computation  of  geoid  undulations  associated  with  wavelengths  shorter 
than  2*  was  done  using  Fourier  techniques  in  a  flat  earth  approximation  as 
described  in  Porsberg  (1985).  The  constant  parts  of  the  topography  and  the 
compensating  masses  were  omitted,  since  they  are  meaningless  in  a  flat  earth 
approximation  (The  geoid  undulation  becomes  infinite).  Then  a  geoid 
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undulation  outside  the  masses  may  be  expressed  by 


Nj  =  (4.2) 

where  y  is  the  normal  gravity,  is  the  potential  from  the  terrain 
d 

Tt  =  GApt  \  \  \  ("*-3) 

A  do 

and  Tc  is  the  potential  from  the  compensating  masses 
d' 

Tc  =  GApt  11^  dzdA  (4.4) 

A  d'o 

G  is  gravitational  constant,  do  is  the  mean  depth  in  an  area  A,  and 
d 

A  first-order  expansion  of  1/r  around  d®  in  eq.  (4.3) 

l/r(2)  =  l/ro-do/ro’(z-do),  to’  =  (xg-Xp)*  +  (yg-yp)*  +  do* 
results  in 

d 

Tj  =  GApt  II  -  ^  (z-do))  dzdA 

A  do 

=  GApt  I  ((^  +  ^)(d-do)  -  (d»-do*))  dA 

=  GApt  I  (d-do)  -  ^  (d-do)*)  dA  (4.5) 

A 

and  for  eq.  (4.4) 

Tc  =  GApc  I  (^^(d'-d'o)  -  (d'-d'o)’)  dA  (4.6) 

A 

The  expressions  for  the  potential  (eq.  (4. 5-4. 6))  are  now  represented  as 
two  dimensional  convolutions,  which  are  suitable  for  the  use  of  Fourier 
techniques,  since  a  convolution  becomes  a  simple  product  in  the  frequency 
domain.  The  series  expansion  of  1/r  that  was  used  in  order  to  obtain  this 
representation,  was  carried  out  to  first  order.  A  zero-order  expansion 
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corresponds  to  a  condensation  of  the  masses  on  a  layer  at  mean  depth.  In  a 
first-order  expansion  more  of  the  geometry  of  the  topography  is  taken  into 
account.  General  expresssions  for  calculation  of  potential  anomalies  Fourier 
techniques  may  be  found  in  Parker  (1972). 

The  accuracy  of  the  expressions  for  the  potential  depends  on  the 
magnitude  of  high-order  terms  that  are  omitted.  According  to  Parker  (1972) 
the  magnitude  of  the  high-order  terms  will  be  smallest,  when  the  reference 
level  is  chosen  so  (d-do)aiax  -  ~(cl~cio)min*  Assuming  that  this  level  is  mean 
level,  this  have  been  done  in  eq.  (4.5)  and  eq.  (4.6).  Furthermore  the 
accuracy  is  highly  depending  on  the  roughness  of  the  topography.  Tests 
carried  out  by  Tziavos  et  al.  (1988)  in  an  area  with  a  very  rough  terrain  show 
that  the  first  and  the  second  term  in  the  expansion  are  sufficient  to  compute 
airborne  gravity  terrain  effects  with  an  accuracy  of  0.25  mgal.  Since  geoid 
undulations  are  smoother  than  gravity  anomalies,  eq.  (4.5)  and  eq.  (4.6)  are 
assumed  to  be  sufficiently  accurate. 

Then  geoid  undulations,  eq.  (4.2),  are  calculated  using  eq.  (4.5)  and  eq. 
(4.6),  and  the  relation:  (d'-d'o)  ~  -o((d-do) 


''l  =  I  {7;  ^  (d-d«)*]  dA 

A 

*  ^  1  (?:  «*-''•>  -  <-*>'  ^  " 

A 

" ;  1  [ik  ■  ■  2  “  Tv) 


(4.7) 


The  integration  in  eq.  (4.7)  are  in  the  form  of  two  convolutions  involving 
(d-do)  and  (d-do)^.  These  convolutions  may  be  performed  most  efficiently  by 
Fourier  transforming  complex  data: 


38 


(4.8) 


+  i 


(d-do)»] 


and  a  complex  kerael: 


[fe  -  -  f^ll 


(4.9) 


Prom  the  computed  spectra  the  apcx^tra  of  each  component  may  be  isolated 
through  the  conjugate  symmetries  of  the  Fourier  transform.  The  kernel  was 
set  to  zero  at  distances  larger  than  2*»  so  only  inner  zone  effects  were 
computed.  If  a  border  zone  larger  than  2*  is  used,  problems  with 
non-periodic  wavelengths  are  avoided  and  windowing  are  not  needed.  A 
maximal  distance  of  2*  was  assumed  to  be  sufficient,  since  only  wavelengths 
shorter  than  2*  were  of  interest.  From  the  Fourier  transforms  of  the 
topography  and  the  kernel,  the  Fourier  transform  of  the  geoid  undulations  is 
obtained  by  a  multiplication. 

In  order  to  obtain  geoid  undulations  associated  with  wavelengths  shorter 
than  2*,  a  high-pass  filtering  was  needed.  This  was  done  in  the  frequency 
domain  before  the  inverse  Fourier  transformation.  As  high-pass  filter,  S(o),  a 
step  function,  Sg  (u),  was  used,  which  has  been  soioothed  by  a  Hanning  taper 
Vi{u)  in  order  to  take  the  frequency  dispersion  into  account  (Mesko,  1984). 
That  is 

So(«)  =  (J 

and 

«(»)  =  (i:  (i  *:) 

lo 

then 

S(m)  =  I  So(u')  W(m'-m)  do' 


for  o<oo 
otherwise 

for  lul  <  Aw 
otherwise 
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0 


for  u^ioo~^u) 

11^^  h  ("^1  (««-4«)<«<(«o+Ao) 

1  otherwise  (4.10) 

where  «•  is  the  frequency)  A«i  is  the  frequency  spacing,  and  <4q  is  the  cut-off 
frequency.  The  cut-off  frequency  is  0.5  cycle/degree#  which  roughly 
corresponds  to  a  harmonic  degree  of  180. 

Then  residual  geoid  undulations,  ANi,  from  the  isostatic  earth  model  were 
obtained  using  an  inverse  Fourier'  transform  of  the  Fourier  transform  of  the 
geoid  undulations  multiplied  by  the  high-pass  filter,  eq.  (4.10). 

Residual  geoid  undulations  were  computed  in  5'  x  5'  "SYNBAPS"  mean 
ocean  depths  was  used.  This  was  done  using  two  compensation  depths: 
De=-50  km  and  0c=-70  km.  If  the  true  compensation  depth  is  assumed  to  be 
30  km,  this  corresponds  to  a  smoothing  of  the  crust/mantle  topography  by 
upward  continuations  of  20  km  and  40  km  respectively.  At  the  locations  of 
the  altimeter  data  the  residual  geoid  undulations  were  found  by  bilinear 
interpolation  in  the  grids. 

4.2  Residual  Terrain  Reduction  of  the  Altimeter  Data  and  the  Local  Empirical 
Covariance  Functions. 

The  residual  terrain  effects  were  f^ompared  with  the  high-pass  filtered 
altimetry  by  evaluating  RMS  values  and  power  spectra  of  the  quantities. 
Furthermore  correlation  coefficients  between  the  two  quantities  were 
calculated.  The  reduction  of  the  altimetry  was  evaluated  in  a  similar  way. 
Finally  local  empirical  covariance  functions  were  estimated  from  the  terrain 
reduced  altimeter  data. 

Correlation  coefficients  between  the  terrain  reduced  altimetry  and  the 
bathymetric  geoid  may  be  used  in  an  evaluation  of  the  earth  model.  They 
cannot  tell  if  a  model  is  correct,  but  they  can  tell  if  a  model  can  be  better.  A 
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zero  correlation  is  obtained,  when  that  part  of  the  observations  that  are 
correlated  with  the  terrain  effect,  are  removed  compleUdy.  If  the  correlation 
coefficient  is  different  from  zero,  the  terrain  model  may  be  used  more 
efficiently.  A  correlation  coefficient  different  from  zero  may  occur  if  wrons 
density  contrasts  are  used.  Then  either  too  leas  or  loo  much  are  removed 
from  the  observations,  which  result  in  either  a  positive  or  a  nef(ative 
correlation  coefficient.  The  use  of  a  wrong  compensation  depth  may  have  the 
same  effect,  since  the  strength  of  the  signal  from  the  compensating  masses 
hereby  changes. 

The  computed  RMS  values  of  the  residual  bethymetric  geoid  undulations 
are  generally  smaller  than  the  RMS  values  of  the  altimetry  (see  Table  4>. 
Furthermore  the  values  associated  with  a  compensation  depth  of  50  km  arcs 
smaller  than  with  a  compensation  depth  of  70  km.  The  computed  correlation 
coefficients  between  the  altimetry  and  the  bathymetric  geoid  (Table  5 1  .show 
that  the  correlations  are  largest  in  Area  3,  where  the  gravity  field  is  rough, 
and  smallest  in  Area  2,  where  the  gravity  field  is  smooth.  A  reason  for 
decreasing  correlations  may  be  that  both  the  altimetry  and  the  ’’SYNBAPS” 
bathymetry  contain  errors.  In  smooth  areas  those  errors  may  dominate  the 
true  signals,  which  results  in  very  small  correlation  coefficients.  Furthermore 
the  isostatic  earth  models  are  highly  idealized  and  many  geological  structures 
are  not  taken  into  account. 


Area 

ANt 

Ah 

Ah-AN,  1 

A 

B 

A 

B 

1 

0.33  m 

0.37  m 

0.53  m 

0.41  m 

0.41  ID 

2 

0.07  - 

0.08  - 

0.45  - 

0.45  - 

0.45  - 

3 

0.86  - 

0.96  - 

1.12  - 

0.42  - 

0.41  - 

Table  4.  RMS  values  of  residual  bathymetric  geoid  undulations,  ANx,  with 
compensation  depths:  Dc=-50  km  (A)  and  Dc=-70  km  (B),  high-pass  filtered 
altimeter  data,  Ah,  and  residual  terrain  reduced  altimeter  data  in  Area  1, 
Area  2,  and  Area  3. 
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Table  5.  Correlation  coefficients  between  altimetry,  Ah,  and  bathymetric 
geuid,  ANi,  with  D^s-SO  km  (A),  and  between  terrain  reduced  altimetry  and 
bathymetric  geoid  with  Dc=-70  km  (B)  in  Area  1,  Area  2,  and  Area  3. 


The  RMS  values  of  the  residual  terrain  reduced  altimetry  are 
approximately  the  same,  and  only  a  slight  change  in  Area  3  is  found,  when  the 
depth  of  compensation  is  changed  from  50  km  to  70  km.  In  Area  1  and  Area  2 
there  are  practically  no  correlation  between  the  terrain  reduced  altimetry  .and 
the  bathymetric  geoid  associated  with  both  compensation  depths.  The 
altimetry  have  therefore  been  reduced  as  well  as  possible  with  both  earth 
models.  In  Area  3  the  terrain  reduced  altimetry  is  still  correlated  with  the 
bathymetric  geoid,  when  a  compensation  depth  of  50  km  is  used.  This 
correlation  is  decreased,  when  a  compensation  depth  of  70  km  is  used.  Sincie; 
the  effects  from  the  masses  above  sea  level  of  Bermuda  are  not  taken  into 
account  when  "SYNBAPS"  data  are  used,  no  further  attempts  to  adjust  the 
earth  model  in  Area  3  were  made. 

The  local  empirical  covariance  functions  and  power  spectra  were  estimated 
from  the  high-pass  filtered  terrain  reduced  altimeter  data,  where  the  residual 
terrain  effects  obtained  with  a  compensation  depth  of  50  km  were  used  in  Area 
I  and  Area  2,  and  the  residual  terrain  effects  obtained  with  a  compensation 
depth  of  70  km  were  used  in  Area  3.  These  results  together  with  the  results 
obtained  before  the  terrain  effects  were  subtracted,  are  shown  in  Figure  7-9. 

The  power  spectra  show  that  the  reductions  are  larger  than  50%  (=3  dB) 
at  frequencies  between  0.75  and  1.75  cycles/degree  in  Area  1,  at  0.75 
cycles/degree  in  Area  2,  and  between  0.50  and  1.75  cycles/degree  in  Area  3 
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Figure  7.  Upper:  Local  empirical  covariance  functions  calculated 
in  Area  1  from  high-pass  filtered  altimeter  data  (o)  and  from 
residual  terrain  reduced  altimeter  data  (+)  (units  in  m* ) .  Lower 
Power  spectra  calculated  in  Area  1  from  lO'xlO'  mean  values  of 
high-pass  filtered  altimeter  data  (o),  residual  terrain  reduced 
altimeter  data  (x),  and  bathymetric  geoid  (+).  The  bathymetric 
geoid  was  calculated  using  a  compensation  depth  of  50  km.  Note 
that  the  symbols  have  changed  in  the  lower  figure. 
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Figure  8.  Upper:  Local  empirical  covariance  functions  calculated 
in  Area  2  from  high-pass  filtered  altimeter  data  (o)  and  from 
residual  terrain  reduced  altimeter  data  ( ■*■ )  (units  in  m*  )  .  Lower: 
Power  spectra  calculated  in  Area  2  from  lO'xlO'  mean  values  of 
high-pass  filtered  altimeter  data  (o),  residual  terrain  reduced 
altimeter  data  (x),  and  bathymetric  geoid  (+).  The  bathymetric 
geoid  was  calculated  using  a  compensation  depth  of  50  km. 
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Figure  9.  Upper:  Local  empirical  covariance  functions  calculated 
in  Area  3  from  high-pass  filtered  altimeter  data  (o>  and  from 
residual  terrain  reduced  altimeter  data  (*)  (units  in  mM .  Lower: 
Power  spectra  calculated  in  Area  3  from  lO’xlO'  mean  values  of 
high-pass  filtered  altimeter  data  (o),  residual  terrain  reduced 
altimeter  data  (x),  and  bathymetric  geoid  (+).  The  bathymetric 
geoid  was  calculated  using  a  compensation  depth  of  70  km. 
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At  higher  frequencies  the  power  of  the  terrain  effects  is  very  small  compared 
with  the  power  of  the  altimetry. 

The  estimated  terrain  reduced  geoid  height  variances  in  the  three  local 
areas  were  0.096  m’,  0.108  m’,  and  0.147  m*  respectively.  The  largest 
reduction  took  place  in  Area  3  where  the  geoid  height  variance  decreased  88% 
(=9.21  dB)t  when  the  terrain  effects  were  subtracted.  In  Area  1  and  Area  2 
the  reduction  were  55%  (=3.47  dB)  and  4%  (=0.18  dB)  respectively. 

The  effects  of  reducing  the  altimeter  data  with  terrain  effects  computed 
from  the  "SYNBAPS"  bathymetry  are  that  the  geoid  height  variances  decreased 
from  0.11-1.22  m^  to  0.10-0.15  m^.  The  magnitude  of  the  unmodelled  parts  of 
the  gravity  field  have  hereby  decreased  and  the  differences  between  the 
geoid  height  variances  in  the  three  local  area  have  decreased  remarkably. 
This  may  be  important  in  a  determination  of  a  covariance  function  model. 

5.  Determination  of  Covariance  Function  Models  from  Local  Empirical  Geoid 

Height  Covariance  Functions 

In  least  squares  collocation  a  covariance  function  model  that  represents 
the  local  empirical  covariance  function,  is  needed.  It  is  therefore  an  important 
step  to  determine  such  a  model.  In  section  2.3  a  method  for  adjusting  a 
Tscherning/Rapp  model  to  empirical  covariance  values  is  described.  In  this 
chapter  the  use  of  this  method  is  described,  when  the  estimated  geoid  height 
autocovariance  values  are  used. 

A  determination  of  a  covariance  function  model  from  a  local  empirical 
geoid  height  autocovariance  function  is  known  to  provide  a  result,  where  the 
depth  to  the  Bjerhammar  sphere  is  not  well  determined  (Knudsen,  1987).  This 
has  been  verified  by  adjusting  a  model  to  the  empirical  covariance  functions 
estimated  from  the  altimeter  data  in  Area  I,  Area  2,  and  Area  3.  Since 
resiaining  long  wavelength  parts  of  the  gravity  field  have  been  removed,  the 
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scale  factor,  a,  was  fixed  to  be  zero.  The  results  showed  that  the  standard 
deviations  of  the  estimated  depths  to  the  Bjerhammar  spheres  were  about  the 
same  magnitude  as  the  parameters  themselves. 

Since  the  depth  to  the  Bjerhammar  sphere  is  not  well  determined  from 
geoid  height  covariance  values  only,  it  was  decided  to  fix  this  parameter  to 
different  depths,  and  adjust  the  factor,  A,  only.  Then  gravity  anomal)' 
variances  were  calculated  from  the  adjusted  covariance  function  models.  The 
results  from  the  adjustments  where  covariance  values  estimated  from  the 
altimeter  data,  before  the  residual  terrain  reduction,  in  Area  1,  Area  2,  and 
Area  3  are  shown  in  Table  6-8  respectively. 


R-H» 

A 

Q 

0.5  km 

1195-10’(m/s)* 

1.02 

0.168  n* 

857  ngal‘ 

1.0  - 

1249 

1.04 

0.166  - 

693  - 

2.0  - 

1359 

1.07 

0.163  - 

545  - 

4.0  - 

1600 

1.13 

0.159  - 

416  - 

8.0  - 

2165 

1.23 

0.152  - 

310  - 

Table  6.  Results  from  the  adjustment  of  a  covariance  function  model  to 
the  empirical  covariance  function  in  Area  1  before  terrain  reduction. 


R-R« 

A 

Q 

CArf. 

0.5  km 

610-10Mni/s)* 

WBM 

0.086  m> 

437  mgal^ 

1.0  - 

639 

0.085  - 

354  - 

698 

0.084  - 

280  - 

827 

0.082  - 

215  - 

8.0  - 

1140 

0.079  - 

162  - 

Table  7.  Results  from  the  adjustment  of  a  covariance  function  model  to 
the  empirical  covariance  function  in  Area  2  before  terrain  reduction. 


R-Hb 

A 

Q 

CAg^  Atf 

8530- 10’ (m/s)* 

0.67 

1.197  m’ 

6111  mgal’ 

1.0  - 

8727 

0.68 

1.162  - 

4842  - 

9506 

0.71 

1.143  - 

3810  - 

11152 

0.77 

1.108  - 

2902  - 

8.0  - 

13463 

0.89 

1.018  - 

2078  - 

Table  8.  Results  from  the  adjustment  of  a  covariance  function  model  to 
the  empirical  covariance  function  in  Area  3  before  terrain  reduction. 

The  results  from  the  adjustments  of  covariance  function  models  to  the 
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empirical  covariance  values  (Table  6-8)  show  that  the  effects  of  a  dramatic 
chanfte  of  the  depth  to  the  Bjerhamntar  sphere  from  0.5  km  to  8.0  km  are 
quite  small,  when  the  Q  values  and  the  model  tfeoid  height  variances  are 
considered.  The  effect  on  the  model  gravity  anomaly  variance,  however,  is 
large.  This  means  that  determinations  of  the  depth  to  the  Bjerhammar  sphere 
may  be  done,  if  estimated  gravity  anomaly  variances  are  available.  Without 
other  information  than  the  geoid  height  covariance  values  the  determination  of 
a  covariance  function  model  is  difficult.  Furthermore  the  gravity  fields  in  the 
three  local  areas  have  different  characteristics,  so  general  assumption  abouth 
the  gravity  fields  cannot  easily  be  made. 

Then  covariance  function  models  were  adjusted  to  fit  the  empirical 
covariance  values  estimated  from  the  residual  terrain  reduced  altimeter  data. 
The  potential  degree  variances  associated  with  the  terrain  reduced  gravity 

field  were  also  assumed  to  tend  to  zero  somewhat  faster  than  i~^,  so  the 

Tscherning/Rapp  model  could  be  used.  It  is  of  course  a  problem  to  decide,  if 
it  is  reasonable  to  model  the  degree  variances  associated  with  the  terrain 
reduced  gravity  field  using  a  Tscherning/Rapp  model.  The  results  obtained 
by  Porsberg  (1986)  in  the  Nordic  countries,  where  accurate  terrain  models  are 
available,  justify  the  use  of  a  Tscherning/Rapp  model.  If  inaccurate  terrain 
models  are  used,  the  gravity  field  will  only  be  reduced  at  wavelengths  longer 
than  the  resolution  of  the  terrain  models,  and  the  shorter  wavelengths  will 

remain  unreduced.  Then  the  covariance  function  model  should  model  the 

degree  variances  in  a  similar  way,  as  if  a  spherical  harmonic  expansion  had 
been  subtracted.  This  may  be  done  by  using  a  sort  of  error  degree  variances 
up  to  a  harmonic  degree  corresponding  to  the  resolution  of  the  terrain  model, 
and  a  Tscherning/Rapp  model  above  this  harmonic  degree. 

The  results  from  the  adjustments  of  covariance  functions  associated  with 
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residual  terrain  reduced  gravity  fields  (Table  9-11)  show  the  same  problems  as 


described  above.  The  imporlatice  of  these  results  are  that  the  differences 
between  the  three  local  areas  have  become  smaller.  The  results  obtained  in 
Area  1  are  very  similar  to  the  results  obtained  in  Area  2.  In  Area  3  the 
respective  variances  are  about  57%  larger  than  the  variances  calculated  in 
Area  2.  This  is  probably  due  to  an  incomplete  modelling  of  the  terrain  above 
sea  level. 

The  major  differences  between  tlie  characteristics  of  the  gravity  fields  in 
the  three  local  areas  have  been  removed,  which  should  make  it  easier  to  use 
general  assumptions  about  the  gravity  field,  when  covariance  func;tion  models 
are  determined.  Furthermore  the  magnitude  of  the  unknown  parts  of  the 
gravity  field  have  decreased,  which  decreases  the  effects  of  using  a  wrong 
covariance  function  model.  In  Area  3  a  change  of  depth  to  the  Bjerhammar 
sphere  from  2.0  km  to  4.0  km  results  in  a  change  in  the  gravity  anoma/y 
variance  corresponding  to  7.9  mgal  (RMS)  without  terrain  correction  and  2.5 
mgal  (RMS)  with  terrain  correction. 


R-Rr  _ 

A 

Q 

CAg- Ag 

0.5  km 

566* 10’ (m/s)* 

0.95 

0.079 

406  mgal* 

1.0  - 

592 

0,96 

0.079  - 

329 

2.0  - 

646 

0.99 

0.078  - 

259 

4.0  - 

763 

1.04 

0.076  - 

199 

8.0  - 

1.16 

0.073  - 

149 

Table  9.  Results  from  the  adjustment  of  a  covariance  function  model  to 
the  empirical  covariance  function  in  Area  1  after  terrain  reduction. 


R-Rr 

A 

Q 

_ C/-./- _ 

CAg. Ag 

0.5  km 

560-10^(m/s)* 

401  mgal^ 

l.O  - 

586 

325 

2.0  - 

639 

256 

4.0  - 

755 

0.075  - 

196 

8.0  - 

1039 

1.29 

148 

Table  10.  Results  from  the  adjustment  of  a  covariance  function  model  to 
the  empirical  covariance  function  in  Area  2  after  terrain  reduction. 
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R-Rb 

A 

Q 

_ _ 

CAtf.Ag 

0.5  km 

881-10=(n/s)* 

0.79 

632  agal> 

i.O  - 

922 

0.81 

512 

2.0  - 

1007 

0.84 

404 

4.0  - 

1192 

0.89 

310 

8.0  - 

1643 

0.97 

234 

Table  11.  Resuils  from  the  adjustment  of  a  covariance  function  model  to 

the  empirical  covariance  function  in  Area  3  after  terrain  reduction. 

Assuming  that  the  decay  of  the  potential  degree  variances  associated  with 
the  terrain  reduced  gravity  field  is  -3.6  corresponds  to  a  covsu’iance  function 
model  with  R-Rb=3.0  km.  In  Area  1  and  Area  2  a  model  with  A=711'10’(m/s)* 
having  a  gravity  anomaly  variance  of  (IS  mgal)*  may  be  used.  In  Area  3  the 
factor  A  and  the  variance  may  be  multiplied  by  1.57. 

In  this  chapter  covariance  function  models  have  been  adjusted  to  the  local 
empirical  covariance  values  estimated  in  the  three  local  areas  from  altimeter 
data.  The  results  show  that  the  characteristics  of  the  terrain  reduced  gravity 
fields,  as  observed  by  the  altimeter  data,  are  quite  similar.  The  determination 
of  covariance  function  models,  however,  still  need  more  information  about  the 
gravity  fields  than  provided  by  the  altimeter  data,  or  generalizations  about 
the  gravity  fields.  In  the  three  local  areas  covariance  function  models  were 
determined  using  the  assumption  that  the  decay  of  the  potential  degree 
variances  associated  with  the  terrain  reduced  gravity  field  is  -3.6. 

6.  Summary  and  Conclusions 

In  this  report,  techniques  for  the  estimation  of  local  empirical  covariance 
functions  are  described.  The  importance  of  removing  those  parts  of  gravity 
field  that  are  associated  with  wavelengths  longer  than  the  extent  of  the  local 
areas,  is  explained  and  a  method  for  doing  it  is  given.  Furthermore  a 
technique  for  the  computation  of  residual  geoid  heights  from  an 
Airy-Heiskanen  isostatic  earth  model  is  described.  Then  the  use  of  satellite 
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altimeter  data  for  the  estimation  of  gravity  field  related  quantities  in  looal 
areas  is  evaluated  and  the  variability  of  the  gravity  field  associated  with 
harmonic  degree  greater  than  180  is  studied  in  three  2*x2*  areas  by 
comparing  local  empirical  covariance  functions  computed  from  altimeter  data 
before  and  after  a  residual  terrain  reduction.  Finally  covariance  function 
models  are  adjusted  to  fit  the  local  empirical  covariance  functions. 

From  the  merged  SEASAT/GEOS-3  altimeter  data  relative  to  0SU81,  geoid 
height  observations  associated  with  harmonic  degrees  greater  than  180  were 
obtained  in  three  steps.  In  the  first  step  the  altimeter  data  were  corrected 
for  time  varying  components  by  a  crossover  adjustment.  In  the  second  step 
remaining  long  wavelength  parts  due  to  errors  in  OSU81,  sea  surface 
topography,  and  correlated  parts  of  remaining  orbit  errors  were  removed  by  a 
highpass  filtering.  In  the  last  step  the  noise  terms  associated  with  the 
altimeter  data  were  modified  in  order  to  obtain  noise  terms  associated  with  the 
altimeter  data  as  geoid  height  observations.  The  results  from  the  crossover 
adjustment  show  that  uncorrelated  parts  of  the  time  varying  components  in 
general  were  well  removed.  Some  short  wavelength  variations  due  by  changes 
in  the  Gulf  Stream,  however,  appeared  not  to  be  well  modelled  in  the 
adjustments.  A  reason  for  that  may  be  that  a  (too)  simple  covariance  function 
associated  with  the  time  varying  components  was  used.  The  covariance 
function  was  originally  designed  to  model  sea  surface  variations  having 
wavelengths  of  about  1000  km  (Knudsen,  1987a),  and  the  variations  due  to 
changes  in  the  Gulf  Stream  have  a  characteristic  wavelength  of  550  km  and  an 
amplitude  up  to  about  25  cm  (Menard,  1983).  A  more  detailed  study  of  the 
non-gravimetric  signals  are  needed,  so  covariance  functions  describing  the 
correlations  in  space  and  time  between  those  signals  can  be  determined. 

The  effects  of  wavelengths  longer  than  the  extent  of  the  local  area  on  an 
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estimation  of  a  local  empirical  covariance  function  were  evaluated  by  comparing 
empirical  covariance  functions  and  power  spectra  estimated  before  and  after 
the  highpass  filtering.  The  effects  were  found  to  be  considerable.  Especially 
the  results  obtained  in  the  "smooth"  Area  2  show  the  importance  of  removing 
those  long  wavelength  parts.  Here  the  magnitude  of  the  long  wavelength 
parts  is  more  than  twice  as  big  as  the  magnitude  of  the  highpass  filtered 
geoid  heights.  In  the  estimation  of  the  long  wavelength  parts  a  smoothing  of 
the  observations  was  done  using  a  rectangular  sine  filter.  A  rectangular  sine 
function  is  a  better  choice  than  a  box-car  function,  but  it  is  not  isotropic.  It 
is  felt  that  the  relations  between  the  spherical  and  the  planar  representation 
have  to  be  studied  more  carefully,  so  more  optimal  low-pass  filters  can  be 
derived. 

Residual  terrain  effects  were  calculated  from  "SYNBAPS"  5'x5'  mean 
bathymetry  using  Fourier  techniques.  As  earth  model  an  isostatic 
Airy-Heiskanen  model  was  used,  but  the  depth  of  compensation  was  increased 
in  order  to  simulate  a  regional  Vening  Meinesz  model.  Residual  terrain 
reduced  altimeter  data  were  computed  using  compensation  depths  of  50  km  and 
70  km  respectively.  The  results  show  that  the  reductions  of  the  altimeter 
data  practically  are  unaffected  by  the  change  in  compensation  depth,  when 
wavelengths  shorter  than  2*  are  considered.  The  important  thing  is  that  the 
RMS  values  of  the  altimetry  have  decreased  to  about  the  same  level.  In  the 
following  tests  a  compensation  depth  of  50  km  was  used  in  Area  1  and  Area  2, 
ans  70  km  was  used  in  Area  3,  since  the  smallest  correlations  between  the 
terrain  effects  and  the  terrain  reduced  altimetry  were  obtained  using  these 
compensation  depths. 

The  results  from  the  comparisons  of  the  local  empirical  covariance 
functions  and  the  power  spectra  calculated  from  the  highpass  filtered  altimeter 
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data  before  and  after  the  residual  terrain  effects  were  subtracted,  show  that 
the  variances'  are  reduced  by  about  3  dB,  0  dB,  and  9  dB  in  three  local  areas. 
As  expected  the  reductions  were  found  to  depend  on  the  roughness  of  the 
gravity  fields,  and  the  estimated  geoid  height  variances  decreased  to  values 
between  0.10  m*  and  0.15  m*,  when  the  terrain  effects  were  removed.  The 
remaining  parts  of  the  gravity  fields,  having  a  RMS  value  of  about  0.35  m,  are 
not  modelled  by  the  Airy-Heiskanen  isostatic  earth  models  that  were  used. 
The  power  spectra  indicate  that  only  a  little  reduction  took  place  at 
frequencies  larger  than  1.75  cyclea/degree.  The  spectra  of  the  altimetry  are 
influenced  by  noise,  but  the  spectra  of  the  bathymetric  geoids  do  decrease 
rapidlly  above  this  frequency.  This  suggests  that  the  resolution  of  the 
"SYNBAPS"  bathymetry  is  about  1.75  cycles/degree  (a  wavelength  of  63  km) 
and  poorer  than  the  resolution  of  the  altimetry,  which  in  chapter  3  was 
determined  to  be  about  3  cycles/degree  (a  wavelength  of  37  km).  Therefore 
the  modelling  of  the  high  frequency  parts  may  be  improved  by  using  a  more 
accurate  bathymetry,  but  changes  in  the  actual  density  distribution  may  also 
cause  high  frequency  changes  in  the  gravity  field. 

The  results  from  the  adjustments  of  covariance  function  models  show  that 
more  information  about  the  gravity  fields  than  provided  by  geoid  height 
observations,  or  strong  generalizations  are  needed,  when  covariance  function 
models  are  determined.  A  removal  of  the  terrain  effects,  however,  removes  the 
major  differences  between  the  gravity  fields,  which  should  make  it  easier  to 
use  general  assumptions  about  the  residual  gravity  field  and,  perhaps  mure 
important,  decrease  the  effects  of  using  a  wrong  covariance  function  model. 
Therefore  the  use  of  terrain  effects  in  gravity  field  modelling  have  reduced 
some  of  the  problems  in  the  determination  of  a  covariance  function  model.  The 
results,  however,  may  indicate  that  the  accuracy  of  the  terrain  model  is  too 
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poor  lo  reduce  the  hi^ch  frequency  parts  of  the  gravity  field,  which  is  not 
taken  into  account  in  the  covariance  function  modelling.  This  explains  why 
the  model  variances  are  smaller  than  the  empirical  variances.  A  study  of  the 
use  of  other  covariance  function  models  may  result  in  a  better  description  of 
degree  variances  associated  with  the  terrain  reduced  gravity  field. 
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